#ifndef IRMM_RAND_H
#define IRMM_RAND_H

#include "Matrix.h"
#include "Vector.h"
#include "assert.h"
#include "mcPublic.h"

#define INF		1e99

void	ztu_srand (unsigned int seed);
int		ztu_rand(void);
unsigned int ztu_seed(void);


class Rand{

public:

	// default constructor generates 1st random number and stores it in r
	Rand();

	// returns a uniform Random deviate between 0.0 and 1.0
	double ran2(int&);
	// Generates a new random number, stores it in r and then returns r
	double NewRand();
	// returns a random number
	double GetRand();
	// acceptance with the specified probability
	inline bool	Accept(double x){if (x<=1e-20) return false; if(NewRand()<=x) return true; else return false;};

protected:

	// seed
	int num;
	// random number
	double r;
};

////////////////////////////////////////////////////////////
// Generates a sequence of random numbers with a particular
// mean value and standard deviation.
class Gauss : public Rand {

public:

/* constructor takes the mean value and the standard deviation of the desired distribution
	then stores the first Gaussian random number in variable g.  Sets
	count to zero. */
	Gauss(const double m = 0, const double s =1);

	double			GetOne(void);							// get one gaussian random variable
	double			GetOne(const double m, const double s);
	mcmcMatrix<double>	GetMatrix(const int m, const int n=1);  // get gaussian random matrix
	/* box_muller function takes two random numbers with uniform distribution
	and returns a random number with Gaussian distribution.
	It also stores the random number in the variable g.the variable g.*/
	double box_muller(const double, const double);

protected:
	double g;
	double mean;
	double sigma;
};


// random field
class RF
{
public:
	// construction & destruction
	RF();
	~RF();

	// operators
	RF&				operator=(const RF &x);
	friend  ofstream & operator>>(RF &x, ofstream &strm);
	friend  ifstream & operator<<(RF &x, ifstream &strm);

	// operatsion
	virtual void	Set(const double x1, const double x2, const double d1);
	virtual void	Set(mcmcVector<float> &votes);
	virtual void	Set(mcmcVector<int> &votes);
	virtual	void	SetOneBin(const double x, const double p);
	virtual	void	AddOneBin(const double x, const double p);

	inline  int	GetCount(void){return count;};
	inline	double	GetDelta(void){return d;};
	inline	double	GetValue(const int index){return x0+d*index;};
	virtual double	pdf(const double x);
			bool	InRange(const double x);
	virtual double	cdf(const double x);
	virtual double	log_pdf(const double x);
	virtual double	pdf_T(const double x, const double T);
	virtual double	d1_log_pdf(const double x);
	virtual bool	Read(const char *lpFile);
	virtual bool	Save(const char *lpFile);
	virtual void	SetCdf(const bool bBig=true);
	inline	double	Sampling(void);
	virtual	double	SamplingNonSetCdf(void);
	virtual int		SamplingIntNonSetCdf(void){return IntClose(SamplingNonSetCdf());};
	virtual bool	normalize(void);
	virtual void	InitValue(const double value);

public:
	mcmcVector<double>	bins;		// data
	int			count;		// count of pdf
	double			x0;			// starting value
	double			d;			// interval

protected:
	__forceinline  int		GetIndex(const double x);

	// data member
	mcmcVector<int>	cdf_bins;
	double			cdf_d;		// cdf interval
private:
	Rand			r_g;
};

__forceinline int RF::GetIndex(const double x)
{
	
	double t1 = (x-x0)/d;
	double y1 = (double)ceil(t1);
	double y2 = (double)floor(t1);
	if (fabs(y1-t1)<fabs(y2-t1))
		return (int)y1;
	else
		return (int)y2;
	/*
	double t1 = (x - x0) / d;
	int y = (int)t1;
	if (t1 > 0){
		if ((double)t1 - y < 0.5){
			return y;
		}else{
			return y + 1;
		}
	}else{
		if ((double)y - t1 < 0.5){
			return y;
		}else{
			return y - 1;
		}
	}*/
	/*
	double t1 = (x - x0) / d;
	int y = ceil(t1);
	if ((double)y - t1 < 0.5){
		return y;
	}
	else{
		return y - 1;
	}
	*/
	/*
	double t1 = (x - x0) / d;
	int y = (int)t1;
	if ((double)t1 - y < 0.5){
		return y;
	}else{
		return y + 1;
	}
	*/
}

inline double RF::Sampling(void)
{
	// sample a uniformly distributed value
	double	v=r_g.NewRand();
	int	index=(int)(v/cdf_d);
	if (cdf_bins.Valid(index))
		return	d*cdf_bins(index)+x0;
	else
		return -1.0;
}

// uniform distribution
class RF_uniform: public RF
{
public:
	RF_uniform();
	RF_uniform(double x1,double x2,double d1);	// u1: mean, w1:width, d1:interval

	virtual void	Set(const double x1,const double x2,const double d1);	// reset pdf function

protected:
	double	u,w;
};

// possion distribution
class RF_gaussian: public RF
{
public:
	RF_gaussian();
	RF_gaussian(double u1,double sigma1,double d1);

	virtual void	Set(const double u1,const double sigma1,const double d1);

protected:
	double u,sigma;
};

// possion distribution
class RF_possion: public RF
{
public:
	RF_possion();
	RF_possion(double lamada);

	virtual void	Set(const double lamda);

protected:
	double lamda;
};

// pose of traffic lights pdf
class RF_pose: public RF
{
public:
	RF_pose();
	RF_pose(double sigma1, double d1);

	virtual void Set(const double sigma1, const double d1);
protected:
	double sigma;
};


///////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////
class RF3D
{
public:
	// construction & destruction
	RF3D();
	~RF3D();

	// operators
	RF3D&	operator=(const RF3D &x);

	// operatsion
	virtual void	Set(const double x1_min, const double x1_max, const double d1_interval,
						const double x2_min, const double x2_max, const double d2_interval,
						const double x3_min, const double x3_max, const double d3_interval);
	inline	void	SetOneBin(const double x1, const double x2, const double x3, const double p);
	inline  void	AddOneBin(const double x1, const double x2, const double x3, const double p);

	virtual double	pdf(const double x1, const double x2, const double x3);
	virtual void	SetCdf();
	inline	void	Sampling(double &x1, double &x2, double &x3);
	inline	void	SamplingNonSetCdf(double &x1, double &x2, double &x3);
	virtual bool	normalize(void);
	inline void		GetValue(double &x1, double &x2, double &x3,
							 const int dim1, const int dim2, const int dim3);

public:
	SimpleMatrix3D<double>	bins;	// data

protected:
	// normalization
	inline  void	GetIndex(int &dim1, int &dim2, int &dim3,
							 const double x1, const double x2, const double x3);

	// data member
	double			x1_0,x2_0,x3_0;			// starting value
	double			d1,d2,d3;				// interval
	mcmcVector<int>	cdf_bins;
	double			cdf_d;					// cdf interval
private:
	Rand			r_g;
};


inline void RF3D::GetIndex(int &dim1, int &dim2, int &dim3, const double x1, const double x2, const double x3)
{
	double t1,y1,y2;

	// the first dimension
	t1 = (x1-x1_0)/d1;
	y1 = (double)ceil(t1);
	y2 = (double)floor(t1);
	if (fabs(y1-t1)<fabs(y2-t1))
		dim1 = (int)y1;
	else
		dim1 = (int)y2;

	// the second dimension
	t1 = (x2-x2_0)/d2;
	y1 = (double)ceil(t1);
	y2 = (double)floor(t1);
	if (fabs(y1-t1)<fabs(y2-t1))
		dim2 = (int)y1;
	else
		dim2 = (int)y2;


	// the third dimension
	t1 = (x3-x3_0)/d3;
	y1 = (double)ceil(t1);
	y2 = (double)floor(t1);
	if (fabs(y1-t1)<fabs(y2-t1))
		dim3 = (int)y1;
	else
		dim3 = (int)y2;
}

inline void	RF3D::GetValue(double &x1, double &x2, double &x3,
							 const int dim1, const int dim2, const int dim3)
{
	x1 = x1_0+d1*dim1;
	x2 = x2_0+d2*dim2;
	x3 = x3_0+d3*dim3;
}

inline void RF3D::Sampling(double &x1, double &x2, double &x3)
{
	// sample a uniformly distributed value
	double	v=r_g.NewRand();
	int	index=(int)(v/cdf_d);
	if (cdf_bins.Valid(index))
	{
		int		dim1,dim2,dim3;
		int	index2;

		index2 = cdf_bins(index);
		bins.GetIndex(dim1,dim2,dim3,index2);

		x1 = d1*dim1+x1_0;
		x2 = d2*dim2+x2_0;
		x3 = d3*dim3+x3_0;
	}
}

inline void	RF3D::SamplingNonSetCdf(double &x1, double &x2, double &x3)
{
	normalize();

	double	v=r_g.NewRand();
	double	sum=0.0;
	int		i,j,k;

	for (i=0; i<bins.Dim1(); i++)
		for (j=0; j<bins.Dim2(); j++)
			for (k=0; k<bins.Dim3(); k++)
			{
				sum += bins(i,j,k);
				if (sum>=v)
				{
					x1 = d1*i+x1_0;
					x2 = d2*j+x2_0;
					x3 = d3*k+x3_0;
					return;
				}
			}
	x1 = d1*i+x1_0;
	x2 = d2*j+x2_0;
	x3 = d3*k+x3_0;
}

inline void RF3D::SetOneBin(const double x1, const double x2, const double x3, const double p)
{
	int dim1,dim2,dim3;

	GetIndex(dim1,dim2,dim3,x1,x2,x3);
	if (bins.Valid(dim1,dim2,dim3))
		bins(dim1,dim2,dim3) = p;
}

inline void RF3D::AddOneBin(const double x1, const double x2, const double x3, const double p)
{
	int dim1,dim2,dim3;

	GetIndex(dim1,dim2,dim3,x1,x2,x3);
	if (bins.Valid(dim1,dim2,dim3))
		bins(dim1,dim2,dim3) = bins(dim1,dim2,dim3)+p;
}


#endif
